History dependence of mechanical properties in granular systems 



ON 

o 
o 

(N 

a 
>—> 

(N 



c 

I 

CO 



c3 



i 

c 

o 
o 



> 
o 
in 



o 

ON 

o 



Shio Inagaki, X 'B Michio Otuski, 2 'Q and Shin-ichi Sasa 3 

'Department of Physics, Kyoto University, Kyoto, 606-8502, Japan 
2 Yukawa Institute for Theoretical Physics, Kyoto, 606-8502, Japan 
:i Department of Pure and Applied Sciences, University of Tokyo, Tokyo, 153-8902, Japan 

(Dated: January 27, 2009) 

We study the history dependence of the mechanical properties of granular media by numerical 
simulations. We perform a compaction of frictional disk packings in a two-dimensional system by 
controlling the area of the domain with various strain rates. We then find the strain rate dependence 
of the critical packing fraction above which the pressure becomes finite. The observed behavior 
makes a contrast with the well-studied jamming transitions for frictionless disk packings. We also 
observe that the elastic constants of the disk packings depend on the strain rate logarithmically. 
This result provides a experimental test for the history dependence of granular systems. 

PACS numbers: 45.70.-n; 46.25.-y; 83.10.Rs; 62.20.D- 



Granular materials have attracted the attentions of sci- 
entists for a long time Despite the extensive efforts, 
our understanding has not yet been deepened enough to 
construct a fundamental law. One major difficulty arises 
from the process dependence of the macroscopic behav- 
iors. It is rather surprising to know that the mechanical 
properties of even static granular media cannot be de- 
termined from all the information of the constituents. 
An example is given by a so-called stress dip, which is 
a local minimum of a vertical pressure distribution just 
under the apex of a sandpile. This stress dip appears 
by local deposition of grains with a point source whereas 
it does not appear by homogeneous deposition with an 
extended source Such process dependence is often 
called as 'memory' which sounds mysterious. A sandpile 
is merely a collection of sand grains between which the 
interaction is rather simple; they interact each other only 
when they are in contact. What attracts us is that only 
granules store information on how they are piled up. See 
Refs. 0, 0, [H as other related memory effects. 

Granular materials have also been of interest as a jam- 
ming system [1,0]- As an idealized case of jamming tran- 
sitions, the packing problem has been studied intensively 
with frictionless particles [1,0, [l(J U|- The important re- 
sult is that the pressure of disk packings is not zero above 
the random close packing Q. On the other hand, with 
frictional particles, the transition point appears between 
the random loose packing and the random close pack- 
ing [HiGjI. It might be naturally conjectured that the 
packing fraction of the transition depends on the packing 
process. However, such a dependence has not been well 
studied. We aim to get insight into the history depen- 
dence of jamming transitions. Furthermore, we explore 
a possibility to describe the history dependence of me- 
chanical properties of frictional disk packings. 

In this Letter, for the purpose of illuminating the na- 
ture of the history dependence, we propose a frictional 
disk packing in a square box without gravity as the sim- 
plest system. We characterize the packing process only 



by its time scale, which is precisely given by the inverse 
of the strain rate of the shrinkage in the packing process. 
By varying this time scale in a systematic manner, we at- 
tempt to observe the process dependence of mechanical 
properties such as macroscopic elastic constants. 

Here, it should be noted that our numerical experi- 
ments are carried out with identical constituents. In or- 
der to examine history dependence in an explicit way, 
it is of the essence to investigate a characteristic fea- 
ture of granular media with persisting with identical con- 
stituents. This makes a contrast with previous studies in- 
vestigating how macroscopic elastic constants depend on 
the material properties of the particles such as the stiff- 
ness, the friction coefficient, and the shape 3, 15, 
The main discovery in the present work is an existence 
of two time scales which are separated over two digits. 
Interestingly, between the two time scales, macroscopic 
elastic properties are expressed as logarithmic functions 
of the strain rate. 



Model: We consider N frictional circular disks in a 
two-dimensional square box. We assume that the di- 
ameter of disk i, denoted by di, is either d or d/lA 
with an equal ratio, and that its mass nii is given by 
n%i = irdfp/4. The interaction between disks appears 
only when two disks are in contact. Suppose disk i at po- 
sition Ti with velocity v.i and angular velocity uii contacts 
with disk j at rj with Vj and ujj . Using the relative posi- 
tion Tij — r.- L — Tj and the relative velocity Vij — Vi — Vj, 
the normal and tangential parts of the relative velocity 
of contacting points, uf 7 - and v\a , are given as vfj = Va ■ n 



and vjj 



ij't— {diUJi + djU)j)/2, respectively, with the 
normal unit vector n = r^/ \rij\ and the tangential unit 
vector t. Let k n and kt be normal and tangential elastic 
constants, r/ n and % be normal and tangential viscous 
coefficients, and a be the Coulomb friction coefficient for 
sliding friction [17|. Then, the normal force and the 
tangential force acting on disk i from disk j are given 



2 



by 



Ffj = mm(\h 



-k n {di/2 + dj/2 - \rij\) + JJnVf 



^l^hsign^,) 



(1) 
(2) 



with h\j 



huh 



+ %vlj. Here, u\ 3 = f^vjj(t)dt is the 
tangential displacement, where to is the time when the 
disks come into contact. The interaction between a disk 
and the side walls of the box is calculated in the same 
way with the one for the interaction between disks. 

The parameters in our simulations are converted into 
dimensionlcss quantities so that d = 1, k n = 1, and p = 1. 
Then, the parameter values we used are as follows, d = 1, 
fi = 0.4, k t = 0.1, and r) n = ry t = 0.35. Note that the 
normal restitution coefficient is approximately estimated 
as 0.6. We restrict our investigation to the systems with 
N = 1000. For a reference to laboratory experiments, for 
example, in the case of disks with maximum diameter of 
10mm, stiffness of 2.5MPa, thickness of 6mm, and mass 
density of 1.6 x g/cm 3 , the time unit in our simulations 
corresponds to about 3.1 x 10 _4 s. 

The time evolution of the position, velocity and rota- 
tion angle of disk i is described by the equation of motion 
with using the forces described above. In numerical sim- 
ulations, we solve a set of equations of motion with em- 
ploying the Adams method with a time step of 2.2 x 10~ 2 . 
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FIG. 1: Illustration of the packing process, (a) Initial config- 
uration: Disks are placed in a box. Open (green online) and 
filled (blue online) circles represent the particles with diame- 
ters of 1.0 and 1.0/1.4, respectively, (b) Schematic diagram 
of the motion of walls: Width of a box, W(t), is reduced as 
a function of time. The box is shrunken isotropically with 
various strain rates. 



We investigate cases with several values of St with 
SW = 5.0 x 10~ 2 fixed. We confirmed that our claim 
in this Letter is independent of the value of SW if it is 
sufficiently small. We stop moving the walls at certain 
time t = t and we wait until the total kinetic energy 
of disks becomes smaller than 10~ 14 so that we can as- 
sume it to be in a static state. This compaction process is 
characterized by the strain rate measured from the initial 
state: 



W(t = 0) - W(t = t) 
W(t = 0)r ' 



(3) 



Note that since the first contact occurs immediately after 
we start a compaction, we assume the initial configura- 
tion as a reference state to determine the strain rate. 
Below we will regard the strain rate e as the control pa- 
rameter of the problem we consider. 

Evidence of the process dependence: We first study 
the process dependence of the final pressure with the 
packing fraction fixed. In this study, we stop the shrink- 
age of the walls when the packing fraction reaches a pre- 
scribed value. 

For comparison, we consider the case /x = 0.0. As 
shown in the inset of Fig. [5J the results with four differ- 
ent strain rates yield one pressure curve as a function of 
packing fraction. Here, the packing fraction of the sys- 
tem is measured in the domain away from the boundary 
at the distance of 2d in order to remove boundary ef- 
fects. We observe that the pressure starts increasing at 
a critical packing fraction <f> c ~ 0.84, which corresponds 
to the so-called jamming transition point It is a re- 
markable property that the critical packing fraction for 
frictionless particles is determined uniquely irrespective 
of packing processes. 

In contrast, under the existence of the tangential fric- 
tion with the friction coefficient of /i = 0.4, the pressure 
curve strongly depends on the packing process as shown 
in Fig. [2] In particular, as we pack slower, the pressure 
curve approaches the one in the case /i = 0.0 (displayed 
by the triangle symbols). In order to extract a quantita- 
tive relation of the process dependence, we measured the 
critical packing fraction (f> c as a function of e. As shown 
in Fig. [3l our result suggests 



As an initial state at t = 0, disks are scattered in a 
square box with no contact. The width of the box W{t = 
0) is determined so that the initial packing fraction is to 
be 0.5. (See Fig. [TJ(a).) We then carry out the following 
compaction procedure with moving all sides of the box 
simultaneously toward the center of the box. Let SW 
be a small value. The width of the box is reduced with 
SW at intervals of time St. The time dependence of the 
width of the box W is illustrated in Fig. [TJ (b) . At the 
moment the width of the box is reduced with the strain, 
centers of disks are shifted with an affine transformation 
with an equal amount of the strain of the box. 



<t> e (ji = 0A,e)-(f> c (n = 0) 



-a<j, log — (4) 



in the regime e' c < e < e' c with e' c <C e' c , where e' c and e' c 
correspond to the inverse of two cross-over time scales r c 
and t' c , respectively. 

We shall present a phenomenological argument to un- 
derstand Eq. ((4]). Let us fix p of the final state. When 
the shrinkage is faster than the change of contact net- 
work and the force balance is realized in this network, 
the packing fraction in the final state does not seriously 
depend on the strain rate. This provides the first cross- 
over time scale t' c . When the shrinkage is slower than 
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this time scale, the disks have enough time to search for 
a more dense network. The searching time r for a config- 
uration with <f> might be proportional to the number of 
force balance configurations S for a given (f>. Here, due 
to the combinatorial complexity of force balance states, 
£ ~ exp(aiV) in the large N limit with a constant a that 
characterizes the number of local configurations satisfy- 
ing force balance conditions. Since it is reasonable to 
assume that a can be expanded in 4> near the packing 
fraction <f/ c at which the shrinkage rate dependence starts 
increasing, £ depends on <p in an exponential manner. 
The consideration leads us to 



r^exp(A(«£-<)), 



(5) 



where A is a constant. This tendency ceases when the 
packing fraction <f) c is close to 4> c {^ = 0), which provides 
the time scale r c as r c = t' c exp(^4(</) c (/i = 0) — <fr' c )) in the 
case p ~ 0. Since the searching time r is related to 1/e 
linearly, we obtain Eq. J1J). 
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FIG. 2: Pressure p versus packing fraction <f> for the system 
with fi — 0.4 and various values of strain rate e, 1.42 x 10~ 3 
(Cross symbol), 3.56 x 10" 4 (Plus symbol), 1.78 x 10" 4 (As- 
terisk symbol)), and 4.45 x 10 (Square symbol)). For com- 
parison, the result in the case /i = 0.0 and the strain rate 
of 4.45 x 10~ 5 is superposed as the triangle symbols. Inset: 
Results in the case /i = 0.0, with various strain rates. The 
symbols coincide with those in the main frame, but all the 
symbols are on one curve. The average is taken over 20 sam- 
ples. 

Experimental method: In many laboratory experi- 
ments, the packing fraction is not easily measurable. We 
thus propose a method by which the process dependence 
can be detected without observing the packing fraction. 
The key idea here is to focus our attention on mechanical 
properties of the packing. Concretely, we fix the pressure 
p in the packing. That is, we stop the compaction when 
the final pressure becomes an prescribed value. In the 
present study, we consider the case p — 0.025. As is seen 
from Fig. [2j the packing fraction might depend on the 
strain rate of the shrinkage, and therefore it is naturally 
expected that Young's modulus E and Poisson's ration 
v of the packings also depend on the strain rate e. 




FIG. 3: Deviation of the critical packing fraction <f> c (fi = 
0.4) from (j> c (fi = 0.0) as a function of the strain rate of the 
shrinkage, where (f> c (n = 0.0) = 0.842 is assumed. Inset: 
Strain rate independence of critical packing fraction 4> c with 
fi = 0.0. 



These elastic constants of disk packings can be mea- 
sured by performing a bi-axial compression test [l8| with 
the final states of the packing process. We measure in- 
crements of stress, 5o~ij, of the system when the strain, 
Se yy = 0.01 is given while Se x = 0.0. Then, assuming the 
linear elastic theory, we estimate the values of E and v. 
Note that we carefully performed the compression test 
so that the quasi-static condition is satisfied during the 
compression. 

The results are summarized in Fig. |H We find that 
the largest Young's modulus is about 1.54 times larger 
than the smallest one. This fact clearly indicates that 
the elastic properties depend on the packing process de- 
spite the material properties of constituents are identical. 
Furthermore, the shapes of curves in Fig. [4] suggests 



E = E - a E log(e), 
v = is + a l , log(e), 



(6) 
(7) 



in the regime e' c < e < e' c 
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FIG. 4: Young's modulus E and Poisson's ratio v (inset) as 
functions of the strain rate e. Young's modulus E = 1 corre- 
sponds to the dimensionless normal stiffness of the constituent 
itself. 40 samples are taken for each r. 

Here we wish to remark on a possibility of realizing 



4 



laboratory experiments with our idea. For example, for 
the experimental system we consider in the paragraph 
Model, the regime 10~ 6 < e < can be realized by 

reducing the size of a box from 46cm to 38cm with com- 
paction time t cx p satisfying 0.3(s) < t cx p < 30(s). This 
operation is accessible in laboratory experiments. Our 
simple setting will provide us a great advantage to inves- 
tigate history dependence of granular media. We expect 
that the logarithmic dependence of E and v on e will be 
examined by laboratory experiments. 

Conclusion and discussion: To conclude, we have 
found the history dependence of the jamming transition 
points. For this phenomenon, the tangential friction is 
found to be essential. Furthermore, we have present 
a clear evidence for the packing process dependence of 
the elastic properties of disk packings with identical con- 
stituents. 



that Y^ijk) ®)k = ^ n f° r eacil As shown in Fig. [5j the 
extent of crystallization decreases rapidly between the 
two cross-over time scales. Therefore, the understanding 
of the tendency to crystallization in slower operations 
might be a heart of future problems. We will study it 
from several viewpoints. 
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FIG. 5: Bond-orientational order parameter t/>6 as a function 
of the strain rate e. 

The most important question left to be answered is to 
uncover the universal nature of the logarithmic depen- 
dence of E and v on e. As discussed above, the sim- 
ple argument suggests the logarithmic dependence of the 
packing fraction on e under the condition that the pres- 
sure p in the final state is fixed. Then, one naively con- 
jectures that the state with larger packing fraction with 
p fixed corresponds to a state closer to crystal in random 
packings. In order to confirm this conjecture, we mea- 
sure the extent of crystallization of the packings with the 
six-fold order parameter of interbond angle defined by 



V'6 



1 



ond 



1 (jk) 



(8) 



where (jk) represents a pair of disks which are in contact 
with disk i next to each other, Q % - k is the interbond angle 
between bond ij and ik, and -/Vbond is the total number 
of such interbond angles. Here, bond ij means the line 
segment connecting the centers of disk i and disk j. Note 
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